/*
 * fft_ops.hpp
 *
 *  Created on: May 26, 2010
 *      Author: z1083248
 */

#if !defined(FFT_OPS_HPP_)
#define FFT_OPS_HPP_

#include <complex>
#include <cmath>
#include <fftw3.h>

// <Debug>
#include <iostream>
#include <fstream>
#include <string>

using namespace std;
// </Debug>

#include "fft_types.hpp"
#include "fft_constants.hpp"

namespace FftLib
{
namespace Common
{

static inline void SinCos(float32 a, float32 *s, float32 *c)
{
   *s = sinf(a);
   *c = cosf(a);
}

static inline std::complex<float32> SinCos(float32 phi)
{
   float32 c, s;
   SinCos(phi, &s, &c);
   return std::complex<float>(c, s);
}

//
// Simple, yet slow (due to non memory locality)...
// http://en.wikipedia.org/wiki/In-place_matrix_transposition#Square_matrices
static inline void squareMatrixTransposeComplex(std::complex<float32> *data,
      uint32 rowSize)
{
   std::complex<float32> tmp;

   for (uint32 n = 0; n < rowSize - 1; ++n)
   {
      for (uint32 m = n + 1; m < rowSize; ++m)
      {
         tmp = data[n + (m * rowSize)];
         data[n + (m * rowSize)] = data[m + (n * rowSize)];
         data[m + (n * rowSize)] = tmp;
      }
   }
}

//
// Simple, yet slow (due to non memory locality)...
// http://en.wikipedia.org/wiki/In-place_matrix_transposition#Square_matrices
static inline void squareMatrixTranspose(float32 *data, uint32 rowSize)
{
   float32 tmp = 0.0;

   for (uint32 n = 0; n < rowSize - 1; ++n)
   {
      for (uint32 m = n + 1; m < rowSize; ++m)
      {
         tmp = data[n + (m * rowSize)];
         data[n + (m * rowSize)] = data[m + (n * rowSize)];
         data[m + (n * rowSize)] = tmp;
      }
   }
}

static inline void convertDataToComplex(const float32 * const inData,
      std::complex<float32> *outData, uint32 inDataLength)
{
   uint32 complexIndex = 0;

   for (uint32 index = 0; index < inDataLength; index += 2)
   {
      outData[complexIndex++] = std::complex<float32>(inData[index],
            inData[index + 1]);
   }
}

static inline void convertComplexToData(
      const std::complex<float32> * const inData, float32 * outData,
      uint32 outDataLength)
{
   uint32 complexIndex = 0;

   for (uint32 index = 0; index < outDataLength; index += 2)
   {
      outData[index] = inData[complexIndex].real();
      outData[index + 1] = inData[complexIndex].imag();
      ++complexIndex;
   }
}

static inline void convertDataToComplex(const float32 * const inData,
      fftwf_complex *outData, uint32 inDataLength)
{
   uint32 complexIndex = 0;

   for (uint32 index = 0; index < inDataLength; index += 2)
   {
      outData[complexIndex][0] = inData[index];
      outData[complexIndex][1] = inData[index + 1];
      ++complexIndex;
   }
}

static inline void convertComplexToData(const fftwf_complex * const inData,
      float32 * outData, uint32 outDataLength)
{
   uint32 complexIndex = 0;

   for (uint32 index = 0; index < outDataLength; index += 2)
   {
      outData[index] = inData[complexIndex][0];
      outData[index + 1] = inData[complexIndex][1];
      ++complexIndex;
   }
}

static inline void printData(const float32 * const data, uint32 dataLength,
      string filePath)
{
   ofstream myfile;
   myfile.open(filePath.c_str());

   for (uint32 index = 0; index < dataLength; ++index)
   {
      myfile << data[index] << endl;
   }

   myfile.close();
}

static inline void printData(const complex<float32> * const data, uint32 dataLength,
      string filePath)
{
   ofstream myfile;
   myfile.open(filePath.c_str());

   for (uint32 index = 0; index < dataLength; ++index)
   {
      myfile << data[index].real() << endl;
      myfile << data[index].imag() << endl;
   }

   myfile.close();
}

}
}

#endif /* FFT_OPS_HPP_ */
